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We present the results of a large-scale numerical study of the equilibrium three-dimensional 
Edwards-Anderson Ising spin glass with Gaussian disorder. Using parallel tempering (replica ex- 
change) Monte Carlo we measure various static, as well as dynamical quantities, such as the autocor- 
relation times and round-trip times for the parallel tempering Monte Carlo method. The correlation 
between static and dynamic observables for 5000 disorder realizations and up to 1000 spins down to 
temperatures at 20% of the critical temperature is examined. Our results show that autocorrelation 
times are directly correlated with the roughness of the free-energy landscape. 
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I. INTRODUCTION 



Systems in statistical physics with rough free-energy 
landscapes, such as spin glasses, are notoriously difficult 
to equilibrate using Markov chain Monte Carlo meth- 
ods. The basic reason for this difficulty is that extremely 
long times are required to surmount free-energy barriers 
and thus to efficiently and fully explore the minima in 
the free-energy landscape. The parallel tempering algo- 
rithm (also known as "replica exchange Monte Carlo") 
[ll-Q, as well as other multicanonical techniques [1, Q 
have partially overcome this problem. However there is 
little understanding of the relevant time scales for these 
methods. The purpose of this work is to shed light on 
the equilibration time of the parallel tempering algorithm 
and to understand how it is related to equilibrium prop- 
erties such as the roughness of the free-energy landscape. 
To this end, we apply parallel tempering to a large en- 
semble of disorder realizations of the three-dimensional 
Edwards- Anderson (EA) Ising spin glass 0, 0] • Our main 
results are that equilibration times vary significantly from 
one disorder realization to another and that this varia- 
tion is correlated with the structure of the free-energy 
landscape. The fact that equilibration times are broadly 
distributed has been previously noted in a different con- 
text (see, for example, Refs. Q and Q) but the extent 
of variation was not fully appreciated because different 
measures of equilibration time were used in past stud- 
ies. The correlation of dynamic time scales for parallel 
tempering and equilibrium properties has not been pre- 
viously recognized, again because of the use of insensitive 
measures of dynamic time scales. In this study we com- 
pute autocorrelation functions of observables to measure 
equilibration times whereas previous studies have primar- 
ily relied on the round-trip time and related quantities. 
Although our study is focused on parallel tempering and 
the Edwards- Anderson Ising spin-glass model, we believe 
that our qualitative results are likely to be relevant both 
to other disordered systems with rough free-energy land- 



scapes, as well as other multicanonical simulation tech- 
niques. 

Parallel tempering was independently introduced by 
Geyer Marinari and Parisi as well as Hukushima 
and Nemoto Q. Furthermore, an earlier algorithm with 
many of the same essential features was introduced by 
Swendsen and Wang [l^. In parallel tempering many 
copies of the system at different temperatures are sim- 
ulated in parallel. Each copy (henceforth referred to as 
replica) is simulated using a standard, single-temperature 
Monte Carlo scheme such as the Metropolis or heat bath 
algorithms [llj . The set of temperatures is chosen so that 
it includes both high temperatures (where rapid equili- 
bration is possible using the single-temperature Monte 
Carlo method) and low temperatures of interest (where 
equilibration is not feasible using the standard algo- 
rithm). In addition to the single-temperature Monte 
Carlo moves, there are also replica exchange moves in 
which replicas at neighboring temperatures swap their 
temperatures. Replica exchange moves permit repli- 
cas to diffuse in temperature space [l^ . The diffusion 
of replicas in temperature space permits an enormous 
speed-up relative to single-temperature Monte Carlo. In- 
stead of directly surmounting the high free-energy bar- 
riers present at low temperatures, barriers are indirectly 
surmounted via the following process: A replica in one 
free-energy well at the lowest temperature diffuses to the 
highest temperature where it easily moves to another well 
and then diffuses back to the lowest temperature. This 
process is called a "round trip" and it allows parallel tem- 
pering to reduce exponential time scales for surmounting 
free-energy barriers to power law time scales [T3| . 

In the context of simple models with rough free-energy 
landscapes it has recently been shown [l^ that the equi- 
libration time of parallel tempering depends strongly on 
the structure of the free-energy landscape. If there is a 
single dominant minimum in the landscape and no first- 
order transition within the range of temperatures, then 
the equilibration time is short. On the other hand, the 
presence of several nearly degenerate minima or the ex- 
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istence of a first-order transition leads to much longer 
equilibration times. These results in simple model land- 
scapes prompted us to look for similar phenomena in the 
Edwards-Anderson model where different disorder real- 
izations may have quite different free-energy landscapes. 
In this study, rather than measuring the free-energy land- 
scape directly, we used the spin overlap (order parame- 
ter) distribution as a proxy. The overlap distribution is 
closely related to the free-energy landscape but is sub- 
stantially easier to measure in simulations. Based on 
the results from simple model free-energy landscapes, 
we hypothesize that when the overlap distribution has 
a complex structure implying many free-energy minima 
the time scale for parallel tempering will tend to be longer 
than when the overlap distribution is simple. Our results 
confirm this hypothesis. 

This study is part of a larger project whose goal is 
to understand the nature of the low-temperature phase 
of finite-dimensional spin glasses. The low-temperature 
phase of the Edwards- Anderson model is poorly under- 
stood and hotly debated. There are several competing 
theories for the low-temperature behavior of the model 
(T3 - l23j and even thou gh m any large-scale investigations 
have been conducted j2a - [26| . a conclusive theory cor- 
rectly explaining all phenomena has not been agreed 
upon. Our recent simulations [27| point toward a two- 
pure-state picture such as the droplet picture [13, [l^ [28l - 
Isoj or the chaotic pairs picture [l^ [2^ |3l| . One moti- 
vation for this work was to verify that the ensemble of 
disorder realizations used in Ref. [13 was indeed equili- 
brated. Therefore, we have measured both equilibration 
times and overlap distributions for all disorder realiza- 
tions and correlated these quantities. We also demon- 
strate that the conservative use of the equilibration cri- 
terion introduced in Ref. [l^l is sufficient to ensure that 
nearly all disorder samples are equilibrated. 

Our results suggest both a method for improving par- 
allel tempering and a warning when using it in spin-glass 
simulations. On the one hand, we find that many dis- 
order realizations have quite short equilibration times. 
Thus it might be useful to implement an adaptive scheme 
where some disorder realization are simulated for shorter 
times than other realizations. State-of-the-art parallel 
tempering simulations of equilibrium spin glasses require 
huge amounts of CPU time because of the difficulty of 
reaching equilibrium and the need for a large ensemble 
of disorder realizations. Thus, such an adaptive scheme 
has the potential for large savings in CPU time. On 
the other hand, the fact that the overlap distribution (a 
quantity directly related to the controversy surrounding 
the low-temperature phase of spin glasses) is significantly 
correlated with the dynamics of the algorithm serves as 
a warning that one must be extremely careful to ensure 
that essentially all samples are well equilibrated in or- 
der to avoid systematic errors in measuring disordered- 
averaged equilibrium properties. 

The paper is organized as follows: In Sec. |TT] we in- 
troduce the model, the parallel tempering algorithm, as 



well as the measured observables. In Sec. IIIII we show 
our results, followed by conclusions. 

II. MODEL AND METHODS 

A. The Edwards-Anderson Ising spin glass 

We study the three-dimensional Edwards- Anderson 
(EA) Ising spin-glass model, defined by the energy func- 
tion 

= - ^ J,jS,Sj . (1) 

The sum is over the nearest neighbors on a simple cubic 
lattice with periodic boundary conditions and side length 
L. The Ising spins Si take values ±1 and the interactions 
Jij are quenched random couplings chosen from a Gaus- 
sian distribution with zero mean and variance one. 



B. Parallel tempering Monte Carlo 

Parallel tempering, also known as replica exchange 
Monte Carlo, is a powerful tool for simulating systems 
with rough free-energy landscapes [l^, [H, [13 with ap- 
plications across many disciplines. To date, it is the most 
efficient method for simulating spin glasses and other dis- 
ordered systems in more than two dimensions at low tem- 
peratures where dynamics are slow. 

The algorithm works as follows: Several copies (here- 
after also referred to as replicas) of the system with the 
same disorder are simulated in parallel. Each replica is 
simulated at a different temperature using a standard 
Markov chain Monte Carlo (MCMC) technique such as 
the Metropolis or heat bath algorithms [llj, IH, H^l • In 
this work we use the heat bath algorithm because of its 
better performance at low temperature. The set of tem- 
peratures is chosen to span both the low temperatures of 
interest, where equilibration is not feasible using single- 
temperature MCMC methods, and high temperatures 
where equilibration is fast. In addition to the MCMC 
sweeps at each temperature, there are also replica ex- 
change moves between replicas. A proposed replica ex- 
change move involves swapping the temperatures of two 
replicas at neighboring temperatures. A replica exchange 
move is accepted with probability p, where 

p = min[l,e(^-'3')(B-B')], (2) 

Here /3 = 1/T is the inverse temperature of one replica 
and (3' the inverse temperature of the neighboring replica, 
and E and E' are the corresponding energies of the two 
replicas. Equation ^ ensures that the entire Markov 
chain, including both single-temperature MCMC sweeps 
and replica exchange moves satisfies detailed balance and 
converges to equilibrium. Replica exchange allows repli- 
cas to diffuse from the lowest temperature to the highest 
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temperature and back again. These round trips enhance 
mixing and greatly reduce equilibration times for rough 
energy landscapes. 

The parallel tempering algorithm has free parame- 
ters that include the set of temperatures and the ratio 
of single-temperature sweeps to replica-exchange sweeps. 
Optimizing parallel tempering therefore requires an ap- 
propriate choice of these parameters. The choice of the 
number of replicas Nj- and their temperatures involves 
the following trade-off: If there are many closely-spaced 
temperatures, the energy distribution between adjacent 
temperatures overlaps strongly and the probability that 
a proposed replica exchange is accepted, see Eq. ([2]), is 
large. This suggests that the more replicas one uses, the 
better apart from the extra computational work involved 
in carrying out single-temperature sweeps of each replica. 
However, the motion of replicas in temperature space is 
diffusive such that the time scale for a round-trip scales 
approximately as the square of the number of tempera- 
tures. At least for simple model systems, parallel temper- 
ing is optimized if the number of temperatures scales as 
the square root of the number of spins [l^ . The average 
time to complete one round trip is often used to charac- 
terize the performance of parallel tempering. Choosing a 
set of temperatures that minimizes the round trip time is 
one of the ways proposed to optimize parallel tempering 

filial [si. 



C. Simulation parameters 

For the parallel tempering algorithm we use Nt ~ 16 
temperatures between T = 0.2 to 2.0 [131 ■ This set of 
temperatures was chosen heuristically in Refs. [23[ and 
[2^ 1 to perform well for L = 6 and 8. The corresponding 
average acceptance fractions for replica exchange moves 
are plotted in Fig. [TJ It is likely that the simulations 
would be more efficient if more replicas were used near 
and above the critical region (Tc ~ 0.96 where the 
acceptance fractions are small. However, we believe that 
as long as the swapping probability is nonzero, our results 
will not change qualitatively. A bottleneck merely slows 
down the diffusion process but does not prevent it. For 
L = 10, the small acceptance fractions suggest that a 
larger number of temperatures might be more efhcient. 
However, our goal is to understand correlations between 
static and dynamic quantities in parallel tempering for a 
given system size and so we have not sought to optimize 
the algorithm for each system size. We believe that the 
qualitative results obtained using our parallel tempering 
parameters would also hold for other reasonable values 
of the parameters. 

By placing the hottest system temperature at T « 2rcj 
we ensure that equilibration happens quickly, and we use 
one exchange sweep for each heat bath sweep (A^ = 
attempted spin updates). To calculate the spin overlap 
[Eq. ^] we use two copies of the system at each temper- 
ature, therefore for each sample we simulate two sets of 
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FIG. 1: (Color online) Average acceptance fraction as a func- 
tion of temperature T for replica exchange moves. Note that 
Tc ^ 0.96. 



replicas independently from each other. 

As stated before, we employ one replica exchange 
sweep after one heat bath sweep. A heat bath sweep 
corresponds to sequentially attempting to update each 
spin at each replica once using the heat bath algorithm. 
An exchange sweep corresponds to randomly choosing 
Nt — 1 pairs of adjacent replicas and proposing exchanges 
for each. Both types of sweeps together make up one par- 
allel tempering sweep. When no confusion arises we will 
call this unit of time simply a sweep. 

For each system size we simulate approximately 5000 
disorder realizations or samples. The simulations are 
equilibrated for at least 2^* sweeps for L = 6 and 2^'' 
sweeps for L = 8 and 10, then measurements arc done 
for the same number of sweeps, see Table |T] for details. 
A fraction of samples required longer runs to meet the 
equilibration criteria discussed in Sec. HID 21 For sam- 
ples with very short equilibration times we performed 
runs to obtain fine-grained autocorrelation information, 
see Sec. HID 21 In total, the data collection took approx- 
imately 140 CPU years on 12-core AMD Opteron 6174 
CPUs. 



D. Observables 

1. Overlap distributions 

The focus of the paper is to relate static equilibrium 
properties of individual spin-glass samples to the dynam- 
ics of the parallel tempering algorithm acting on that 
sample. The primary quantity that we measure to study 
both equilibrium properties and the dynamics of parallel 
tempering is the spin overlap q. 



q = 



1 ^ 
N ^ 



(3) 



TABLE I: For each system size L we equilibrate and then 
measure for at least 2' Monte Carlo sweeps. Tmin (Tmax) is 
the lowest (highest) temperature used, Nt is the number of 
temperatures, and N^a, is the number of disorder realizations. 
For some L = IQ samples longer runs had to be performed to 
ensure equilibration. 



where N — is the number of spins and the super- 
scripts (1) and (2) indicate two independent copies of 
the system with the same disorder. The thermal and dis- 
order average of the overlap is an order parameter for 
the system and the probability distribution of the ther- 
mally averaged overlap Pj{q) for a given sample J reveals 
aspects of the free-energy landscape for that particular 
sample. The behavior of Pj{q) for large systems and 
at low temperatures is one of the major open questions 
in the theory of spin glasses. Overlap distributions vary 
widely from sample to sample, as can be seen in the three 
examples shown in Fig. [21 Although there is no direct 
mapping between the free-energy landscape and Pj(q), 
it is clear that numerous peaks in Pj{q) imply numerous 
minima in the free-energy landscape. For example, sam- 
ples corresponding to Figs. HJb) and HJc) have a more 
complicated free-energy landscapes than Fig. [2lja). Our 
central hypothesis is that samples with a more compli- 
cated free-energy landscape tend to have longer dynamic 
time scales. 

From Pj{q) we define Ij{qa) for different qo, 



Ijiqo) 



90 



(4) 



<?0 



Ij{qo) is the weight of Pj{q) in an interval near the origin 
and it serves as an approximate measure of the complex- 
ity of Pj{q)- For example, Ij (0.2) « for the samples 
shown in Figs. [^Ja) andjU^b), while /j(0.2) is large for 
the sample shown in Fig. [2][c). We compute /j(qo) for 
eight values of go from go = 0.2 to 0.9. It is worth noting 
that Jj (50 )~especially with go = 0.2-has been extensively 
used in studies of the low-temperature equilibrium prop- 
erties of the EA model [H, [2^. 



2. Characteristic time scales 

We measure two time scales for parallel tempering from 
the autocorrelation function of the overlap at the lowest 
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FIG. 2: Examples of overlap distributions Pj{q) for three 
different disorder realizations J for L = 8. While the distri- 
bution in panel (a) has only two peaks, the distribution in 
panel (c) implies a very complex free-energy landscape. All 
panels have the same horizontal scale. 



temperatures. For an arbitrary observable A, the auto- 
correlation function Tj^{t) is defined via 



{AiO)Ait)) - {Ay 
{A^) - (A)' 



(5) 



Here, A{t) is the observable measured at time t and (■ • • ) 
represents a thermal average [38| . Integrated and ex- 
ponential autocorrelation times are computed from the 
autocorrelation function. The integrated autocorrelation 
time is the integral of the autocorrelation function or, 
for discrete time, the sum of the autocorrelation function. 



E 



(6) 



The integrated autocorrelation time is the time needed 
for two subsequent measurements of A to decorrelate. It 
is a lower bound on the time needed to reach equilibrium 
from an arbitrary initial condition. To calculate t-^^ care 
must be taken when truncating the sum in Eq. ^ to 
avoid large statistical errors (explained below). 

Figure [3] shows two typical examples of the autocorre- 
lation function for the absolute value of the spin overlap 
\q\. For large times, the autocorrelation function is even- 
tually dominated by noise. The noise floor is indicated 
in Figs. [3] by the horizontal (red) solid line. Because the 
noise floor is determined primarily by the number of data 
points used to compute the autocorrelation function, it 
can be chosen to be the same for all simulations of the 
same length. Our truncation procedure is to sum the 
autocorrelation function until it first hits the noise floor. 
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FIG. 3: (Color online) Autocorrelation function for the abso- 
lute value of the order parameter T\q\ {t) as a function of simu- 
lation time for two different disorder realizations and system 



size L 
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represents the noise floor. The diagonal solid (green) lines 
represent exponential fits to the data. Both panels have the 
same horizontal scale. 



We obtain the noise floor for each run length by visual 
inspection of Fi^i [t) for many samples. We note that this 
approach introduces some slight systematic bias in our 
estimates of the integrated autocorrelation times. How- 
ever, the main emphasis of this paper is to understand 
the very large sample-to-sample variations in time scales 
and, compared to these variations, the errors introduced 
by our truncation protocol are small. 

We have computed autocorrelation functions for both 

\q\ and q and from these we have obtained rj''| and r^^, 
respectively. To do this, we first compute a Fourier trans- 
form of the data and then invert the function to ob- 
tain the autocorrelation function. This method is con- 
siderably faster than a direct calculation. For L = 6, 
the simulations are run for 2^^ sweeps, recording data 
for the computation of the autocorrelation functions ev- 
ery 10 sweeps. The noise floor for these runs is set to 
F = 0.005. Some samples have Tj'^| values of order 10 or 
less. For these, shorter data collection runs of 2^* sweeps 
have been done starting from an equilibrated spin con- 



figuration stored at the end of the longer run. For these 
short runs data are collected every sweep to accurately 
measure the autocorrelation function up to time 10 and 
the full autocorrelation function is patched together from 
the short and long runs. The noise floor for these shorter 
runs is set to F = 0.01. If the noise floor is reached be- 
fore t = 10, only the autocorrelation function generated 
by the short run is used. For both L = 8 and L = 10 
the simulations are run for at least 2^'' sweeps, recording 
data every 100 sweeps with shorter runs also needed for 
samples with small Tj'^| values. The noise floor used is 
the same as for L = 6. 

The autocorrelation function is always calculated up 
to a maximum time lag of 1% of the total run length. 
Therefore, for some samples and, in particular for L = 10, 
the noise floor was not reached by this time and longer 
runs were necessary to obtain good statistics. For these 
L = 10 samples we simulate up to 2^^ sweeps. 32 sam- 
ples (approximately 0.6%) stayed above the noise floor 
even for the longest runs. To prevent a biasing of the re- 
sults, these were not included in the analysis. However, 
it would be interesting to study these samples in detail 
in the future. 

There is another method for truncating the sum in 
Eq. ^ introduced by Madras and Sokal [331 where the 
upper limit of the sum is determined recursively. An ini- 
tial upper limit is chosen and is computed for that 
value, the new upper limit is determined as some fac- 
tor (6 is suggested in Ref. [3§]) times the current Tj^^. 
The estimated value of is obtained when this process 
has converged. We experimented with this approach but 
did not flnd that it converged in a consistent way across 
all samples with very different autocorrelation functions. 
However, we note that our approach should have system- 
atic errors no larger than those resulting from the method 
of Madras and Sokal with the upper limit set to 6rj''| be- 
cause the noise floor is not reached until at least IOtj^^I 
except for the very few samples with extremely small Tj'^| 
values. 

The second time scale that we consider is the exponen- 
tial autocorrelation time r^p, defined by the asymptotic 
exponential decay of the autocorrelation function, 

TA{t) - Ae 



-t/T 



(7) 



Except for observables orthogonal to the slowest mode of 
the Markov chain, is the characteristic time of the 
slowest mode and is thus the exponential time scale to 
reach equilibrium from an arbitrary initial state [40| . 

In principle, r^p is the most important time scale 
for studying disordered systems because the main diffi- 
culty is reaching thermal equilibrium (controlled by r^p) 
rather than collecting enough uncorrelated data at equi- 
librium (controlled by r-^t)- Unfortunately, calculating 
T^p for a large number of samples is difficult because it 
requires an automatized fitting process. For some sam- 
ples, such as the one shown in Fig. 3(a) the autocorre- 



lation function decays nearly exponentially over most of 
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FIG. 4: Logarithmic histogram of the integrated autocor- 
relation times Tj'j for for L = 8. The data are normally 
distributed and show that a small fraction of samples is ex- 
tremely difficult to equilibrate. 

the observable time range, allowing for a precise fitting 
and extraction of Tj'^|. However, for many samples the 
autocorrelation function is not exponential for the mea- 
sured times, as can be seen in Fig. |3(b)| and, for these 
samples, an automatic fitting procedure would not be re- 
liable. For that reason, the time scale that we correlate 
with the static properties of the EA spin glass is the in- 
tegrated autocorrelation time. We have studied rixp by 
hand for a small subset of the data set and find that for 
those samples where a reasonable fit is possible, rixp is 

in the range of 1 to 15 times r;''|. Most, but not all, sam- 
ples share the behavior seen in the examples in Figs. [3] 
that there is an initial sharp decline before the (approx- 
imately) exponential decay sets in. This sharp decline 
explains why rixp is typically larger than Tj'jj|. 

Finally, we also measured the round-trip time trt for 
each sample, defined as the average time after equilibra- 
tion is achieved, for a replica to diffuse from the lowest 
to the highest temperature and then back to the lowest 
temperature again [T2j . 



III. RESULTS 

Figure |4] shows the probability distribution for rj^j for 
L = 8. In agreement with previous studies d, [ot 125|. 
the figure reveals that the equilibration times for spin 
glasses are very broadly distributed. The distributions 
for L = 6 and 10 are similar although with different 
means and standard deviations, as shown in Tables HI] 
and mil respectively. Tables HI] and IIIII also show the 
mean and standard deviation of Tj^^j. The values for Tj^^ 

are similar to those for rj'j^|, however the standard devi- 
ations are somewhat smaller. Figure [S] shows a scatter 
plot with each point representing one sample for L = 8. 




FIG. 5: (color online) Scatter plot of t^I^^^ vs. rf^l for L = 8. 
The diagonal (red) solid line corresponds to t^^^ — T^^j.. Note 
that T; is often smaller than Tj^^^ because the decorrelation of 
q requires global spin flips while the decorrelation of \q\ does 
not. 



TABLE II: 
for system 


Mean values of the lo 
sizes L = 6 - 10. 


2:arithms of Tj^^ , 


tI^I and TRT 


L 


6 


8 


10 


loglo(-f]nt) 

logio(^i'.?i) 
logio(-rRT) 


1.7164(58) 
1.6035(90) 
3.6217(28) 


2.8366(75) 
2.8065(94) 
4.4571(34) 


3.8711(84) 
3.8881(97) 
5.2421(42) 


TABLE III: Standard deviations of the logarithms of Tj^^, t-^^I 
and trt for system sizes L — 6 - 10. 


L 


6 


8 


10 


logio(^hit) 
logio(^iLI) 

loglo(TRT) 


0.408 
0.631 
0.194 


0.534 
0.672 
0.242 


0.592 
0.685 
0.293 



The X and y coordinates are determined by Tj^^ and r- '| , 

respectively. The data illustrate that Tj^^ and r-'l are 
strongly correlated, especially for those samples with the 
longest integrated autocorrelation times. The observa- 
tion that Tj'^l can take smaller values than Tj^^ is presum- 
ably a consequence of the fact that the decorrelation of 
q requires global spin flips while the decorrelation of \q\ 
docs not. 

Figures [BJa) and (HKb) show scatter plots of Tj'^| vs 
Ijilo) for L ~ 8 and go = 0.2 and go = 0-8, respectively. 
Each point represents one sample. These figures show 
that Tj'l is correlated with /j(go) and the correlation is 
stronger for larger values of go- H is also striking that 
large values of //(go) arc never associated with very small 
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FIG. 6: Scatter plot for all L = 8 samples of r^^l vs Ij{qo) 
for (a) qo — 0.2 and (b) go ~ 0.8. The empty area underneath 
the scatter for large Ij{qo) clearly shows that realizations with 
large values of /j(go) can never have small values of T;'^. 

integrated autocorrelation times (similar results are ob- 
tained for the other system sizes studied). This suggest 
that when the free-energy landscape is complex, time- 
scales are large, as one would naively expect. 

To quantify the correlations between Tj'J and Ij{qo), 
we use the Pearson correlation coefficient r [4l|. whose 
estimator for a finite data set {{xi,yi)} is 

M 

^{xi - x){yi - y) 
I (8) 

\ 1=0 1=0 

where M is the size of the data set and the overbar indi- 
cates the average over the data set. If r = 1 there is an 
exact linear relationship between x and y while a small 
value of r indicates the absence of any linear dependence. 

Pearson r values for logio(Ti'nt) vs logiQ[/,/(q'o)] for our 
data are shown in Table HVl and Fig.[7]for the system sizes 
studied and different values of go- These results clearly 
demonstrate that longer parallel tempering time scales 
are correlated with more complicated overlap distribu- 
tions. The strength of the correlation is slightly weaker 
for larger systems. The fact that r increases with sug- 
gests that the presence of additional peaks in Pj{q) in- 
creases Tj'^l independent of whether those peaks are near 
the origin or not. 



TABLE IV: Pearson correlation r between the overlap weight 
logjQ[/j(qo)] for various go and the logarithm of the integrated 
autocorrelation time for the overlap logjo(T"int) foi' the three 
system sizes, L. The bottom row shows the correlation coeffi- 
cient between the logarithm of the round trip time logjQ(rRT) 
and logwiAut)- 



qo 


L = 6 


L = 8 


i = 10 


0.2 


0.5299(94) 


0.5355(95) 


0.4889(95) 


0.3 


0.5634(86) 


0.5693(79) 


0.5264(93) 


0.4 


0.5982(80) 


0.5983(75) 


0.5590(91) 


0.5 


0.6286(77) 


0.6217(73) 


0.5920(80) 


0.6 


0.6671(72) 


0.6494(70) 


0.6272(72) 


0.7 


0.7062(63) 


0.6788(66) 


0.6535(65) 


0.8 


0.7385(46) 


0.6975(63) 


0.6752(57) 


0.9 


0.7397(58) 


0.6583(85) 


0.6410(83) 


logio(TRT) 


0.014(19) 


0.114(17) 


0.126(20) 



On the other hand, we find that the round trip time is 
not correlated with tI^I or, by extension, Ij{qQ). Figure 

|8]is a scatter plot of trt vs Tj'^I for L = 8 (compare with 
Fig. [S]). Each point represents one sample and the diag- 
onal line is trt ~ t-^^^I- It is clear from this figure that 

trt has a much smaller variance than Tj'^I and that the 
correlation between the two is minimal. The last row of 
Table HVl shows the Pearson r value between logiQ(TRT) 
and logio('''i',it|) ^^'^ quantifies the lack of correlation be- 
tween trt and tI^^. Note that the authors of Ref. [2^ 
consider a time scale related to the motion of replicas 
in parallel tempering and also find no correlation with 
Ij{qo)- 

If the highest temperature is chosen such that its mix- 
ing time is one heat bath sweep, then trt > '''int because 
there will be no memory of the spin state of the replica 
between successive visits to the lowest temperature, as 
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FIG. 8: (Color online) Scatter plot of the average round 
trip time trt vs t^^I for L — 8. The diagonal (red) line is 
trt = tI^I- The data show almost no correlation between 
these quantities. 
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FIG. 9: (color online) Average link overlap [{qi)]a.v (red cir- 
cles) and 1 — T\u\/d (see Eq. (Illf) . blue squares) as a function 
of time t for L = 8 and T = 0.20. For t > lO'' sweeps both 
data sets agree suggesting that the system is in thermal equi- 
librium. The vertical (green) line in the figure represents the 
time for which 99% of the samples have values of rj''^ less than 
this time. Error bars are smaller than the symbols. 



can be seen in Tables lU and IIIIl as well as Fig. [51 Figure 
|8] shows that trt is much larger than rj'^| for most sam- 
ples. One explanation for this is that (for most samples) 
the spin state of the coldest replica can be decorrelated 
by diffusion from the lowest temperature to an intermedi- 
ate temperature and back; full round trips are not needed 
for decorrelation. 

Note that if we repeat the study for Tj^^ instead of Tj^^I , 
the correlation values would show the same qualitative 
trend as for t-J^I, but they would on average be slightly 
smaller. 

Finally, we also test equilibration with the method in- 
troduced by Katzgraber et al. in Ref. [l^l for short-range 
spin glasses with Gaussian disorder between the spins. 
In this method, the average energy per spin 

-^^[Jij{siSj)]^^, (9) 

where [• • • ]av represents a disorder average and, as before, 
(• • • ) denotes the Monte Carlo average for a given set of 
bonds, can be related to the link overlap 

via an integration by parts over the interactions Jij be- 
tween the spins. In Eq. Nt = dN with d = 3 the 
space dimension represents the number of bonds in the 
system. One obtains: 

r, T\U\ 

[(9;>]av = l-^, (11) 
where T is the temperature. 



Figure [5] shows representative results for L = 8 and 
T = 0.2, the lowest temperature simulated. At approxi- 
mately 10^ Monte Carlo sweeps both data computed di- 
rectly from the energy per spin (blue squares) and data 
computed from the link overlap (red circles) agree within 

error bars. Note that 99% of the t-^I values for this sys- 
tem size are smaller than the point marked with the solid 
vertical green line. Thus, for the simulated time of over 
10^ Monte Carlos sweeps, the data are well equilibrated. 
This suggests that a conservative use of the equilibration 
criterion developed in Ref. [l^ will guarantee that almost 
all of the samples are equilibrated. 

IV. CONCLUSIONS 

We have found one factor that explains the broad dis- 
tribution of time scales in parallel tempering for spin 
glasses: The roughness of the free-energy landscape of 
each individual sample directly affects its equilibration 
time scales. In accordance with previous observations 
on the parallel tempering method [13| , if there are many 
minima in the free-energy landscape, a sample requires 
more time to equilibrate. The equilibrium distribution of 
the spin overlap serves as proxy for the complexity of the 
free-energy landscape. Our results show that individual 
samples need to be tested individually to ensure proper 
equilibration, as previously suggested in Ref. @. How- 
ever, we show that the equilibration test developed in 
Ref. [2^ is a viable alternative for system sizes currently 
accessible in simulations if applied conservatively. Fur- 
thermore, the parallel tempering round-trip times seem 
to be unaffected by the complexity of the free-energy 
landscape and should therefore not be used as a mea- 
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sure of how well a system is equilibrated. 

The clear correlations between algorithmic (nonequi- 
librium) time scales and the equilibrium complexity of 
the free-energy landscape opens the door for alternate 
studies of the nature of the spin-glass state, a problem 
that remains controversial. 
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